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Abstract 

We present an algorithm for Monte Carlo simulations which is able to 
overcome the suppression of transitions between the phases in compact U(l) 
lattice gauge theory in 4 dimensions. 



1. Introduction 

Monte Carlo simulations of compact U(l) lattice gauge theory in 4 dimensions 
near its phase transition are hampered by the strong suppression of transitions 
between the phases. On larger lattices conventional algorithms are not able to 
induce transitions at all. 

We obtain a more appropriate algorithm by using the Wilson action supplemented 
by a monopole term [1] 

S = ^ J2 (l-cos0^.,.) + A^|M,,.| 
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where Mp^^ is the monopole content of 3D cubes [2]. The strength of the first order 
transition decreases with A, the transition ultimately getting of second order [3]. 
Therefore, by making A a dynamical variable, we no longer rely on the tunneling 
between the phases but get a much easier pathway along the mountains of the joint 
probability distribution P(£',A), where E denotes the plaquette energy (Figure 2 
gives an example of P{E^ A)). Our algorithm has the further virtue to be vectorizable 
and parallelizable. 

Before running the dynamical-parameter algorithm one has to determine the phase 
transition line (shown in Figure 1) and the additional term in the action related 
to the prescribed distribution of A. It is crucial that this can be done with low 
computational cost also on larger lattices. Here we present a method by which this 
is achieved. 

2. Outline of method 

Conventional methods simulate the probability distribution 

/iA(0) = exp{-Sx{e))/Z, 

where A is a fixed parameter. In order to make A a dynamical variable we consider 
fi\{Q) as the conditioned probability to get a configuration given a value of A and 
prescribe a probability distribution /(A) to get the joint probability distribution 
/i(0. A) = f[X)fi\[Q). To simulate /i(0. A) we need it in the form 

/i(0,A) = exp(-5'(0,A))/Z (1) 

where 5'(0, A) = 5'a(0) + (/(A) and g{X) is determined by /(A). 

In our application of the algorithm each update of the link variables 0^,^; is followed 
by an update of A, which is done in a discrete set of values Xq with (j' = 1, . . . , n. The 
individual update steps are Metropolis steps in both cases. For the A update the 
proposal matrix ^{Sq^i^qf-\-Sq^qf^i-\-Sq^iSqf^i-\-Sq^nSqf^n) ^ud thc acccptaucc probabilities 
min(l, exp(S'(0, Xq) — 5'(0, Ag/))) are used. 

For each value of q one needs /3(Ag) and g{Xq). We require /3(Ag) /3u,(Ag) where 
is the /3 value where both phases are equally probable. Our condition for fixing 
g{Xq) is /(A) saconst . In order to determine the sets of /3(Ag) and g{Xq) we use the 
fact that in a simulation the transition probabilities between neighboring values of 
A are very sensitive to these quantities. 
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3. Transition probabilities 



To derive relations which can be used in the envisaged determination of /3(Ag) 
and g{Xq) we use the probabihty for the transition from a value Xq to a neighboring 
value Xqi 

W{Q, q- q') = i min(l, exp(^(0, Xq) - ^(0, A,,)) (2) 
and note that detailed balance implies 

/(A,_i)/iA,_, (0)W/(0, q-l;q) = /(AJ/iA,(0)W/(0, 9; 9 - 1) • (3) 
Introducing the average transition probability for a set K{q) of configurations 

PK{q: q') = E /^A,(0)W/(0, q- q') , (4) 

WK{q) 

where wk is the weight of this set, we get 

f{Xq_i) WK'iq - 1) PR'iq -l;q) = f{Xq) WK{q) PR'iq; q - I) ■ (5) 
from averaging (3). 

We now apply (5) to sets of configurations in the hot phase and in the cold phase 
separately. Because we are interested in cases where transitions between the phases 
are extremely rare, sets of this type with numbers of configurations sufficient for the 
present purpose are easily obtained. The respective equations can be considered to 
be independent. Using (5) for K = Kc and K = we get a pair of equations 
which simplifies to 

PKrXq - 1;?) = PKrXq;q- i) 

PKn{q - 1; = PKn{q; 9-1) (6) 

ii fi = fiw and /(A) =const. This is what we exploit to determine fi{Xq) and g{Xq). 

Our strategy is to adjust fi{Xq) and g{Xq) for known /3(Ag_i) and (/(Ag_i) such that 
(6) holds. Starting from given /3(Ai) and arbitrarily chosen g{Xi) in this way we can 
obtain /3(Ag) and g{Xq) for q = 2^ . . . ^n. 

4. Determination of (3{\) and ^(A) 

To start our procedure we use a value of Ai in the region where the peaks of the 
probability distribution related to the phases strongly overlap so that tunneling is 
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no problem and /3{Xi) can easily be obtained by a conventional simulation. Because 
only the differences (/(Ag_i) — g{Xq) are relevant we can choose g{Xi) = 0. Then 
for (j' = 2, . . . , n we consecutively determine /3(Ag) and g{Xq) for known /3(Ag_i) and 

Within a step from q — 1 to q we first get a rough approximation of /3(Ag) by 
extrapolation from former values and an approximate new Xq in about the same 
distance as in former steps. Then we generate the sets of configurations Kc and 
Kh by short Monte Carlo runs with cold and hot start, respectively. Because for a 
set K{q) the quantity 

= ttV^ E W{e,q;q'), (7) 

where Nk is the number of configurations, approximates (4), this allows us to cal- 
culate approximate quantities which should satisfy (6). These are two equations 
determining the two unknown values /3(Ag) and g{Xq). We obtain good estimates 
for them though only approximate quantities enter (6) because the peaks related 
to the phases vary only little with /3. In addition, p[q — l;q) is used to adjust the 
distances between neighboring A values such that one has roughly equal transition 
probabilities for all steps. 

After a larger number of q steps the errors may accumulate. Therefore, we perform 
short runs of the dynamical-parameter algorithm to test wether it does indeed travel 
along the mountains of the distribution in the hot as well as in the cold phase. If 
it gets stuck we slightly increase or decrease the /3(A) in the region of A where 
the transitions fail. We then determine the corresponding values of g{Xq) from the 
conditions 

PKcii - 1; + PR'hii - 1; = PKcii; 9 - i) + pkM', 9-1) (8) 

and try again to run the dynamical algorithm. Typically one or two trials are 
sufficient. 

After performing the actual simulations with dynamical A, improved /3(Ag) can be 
obtained by reweighting [5] the distribution at A values where deviations from equal 
weights of the phases occur. The related new g{Xq) then are obtained from (8). The 
values g{Xq) can be improved by replacing them by g{Xq) -\- ln(/(Ag). 

5. Numerical results 

Figure 1 shows the location of the maximum of the specific heat jSc as function 
of A. The method presented has enabled us to determine this phase transition line 
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Figure 1: Location of phase transition as function of A for 8"* (circles) and 16'^ 
(crosses) lattices. 

also for lattice size 16'^ needing only a small amount of computer time (which by 
conventional simulations due to the suppression of transitions is not possible at all). 
To get the sets of configurations Kc and we used about 2 X 10"* sweeps per A and 
in the short test runs with the dynamical algorithm 4 X 10"* sweeps in total. The 
calculations have been done for n = 25 values of Xq ranging from A = to A = 0.4. 

In our simulations with dynamical A, in which we collected about 10^ sweeps, it has 
become possible to have transitions between the phases also on the 16'^ lattice. They 
confirmed the values which we have found for /3u,(Ag) by the procedure described 
above. Figure 2 shows the (reweighted) distribution P{E^ A) which we obtained 
in our simulations (by the dynamical-parameter algorithm except for A = 0.5 and 
A = 0.6). 

On the 8"* lattice there is still substantial tunneling. This is refiected by the overlap 
of the peaks in the distribution P{E^ A) which we have presented in [3]. There the 
g{Xq) have been determined by a simpler method based on (8). 

The tunneling times between the phases in our algorithm with dynamical A on 
the 8"* lattice are greatly reduced as compared to those of a conventional Metropolis 
algorithm [3]. On the 16'^ lattice, where for the conventional algorithm one observes 
no transitions at all, for A = we get a time of the order of 10"^ for our algorithm. 

For further reduction of the autocorrelation times in addition to avoiding tunneling 
one would have to replace the local Metropolis steps for by more efficient ones 



5 




Figure 2: Distribution P{E^ A) at the phase transition hne on 16'^ lattice. 

(corresponding to the cluster steps in the analogous algorithm for the Potts model 
[6]). 
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